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ABSTRACT 

Data interpolation is a fundamental problem in many applied mathematics and scientific computing fields. 
This paper introduces a modified implicit local radial basis function interpolation method for scattered 
data using polyharmonic splines (PS) with a low degree of polynomial basis. This is an improvement to the 
original method proposed in 2015 by Yao et al.. In the original approach, only radial basis functions (RBFs) 
with shape parameters, such as multiquadrics (MQ), inverse multiquadrics (IMQ), Gaussian, and Matern 
RBF are used. The authors claimed that the conditionally positive definite RBFs such as polyharmonic 
splines r°” Inr and r?"*! ” failed to produce acceptable results”. 

In this paper, we verified that when polyharmonic splines together with a polynomial basis is used 
on the interpolation scheme, high-order accuracy and excellent conditioning of the global sparse systems 
are gained without selecting a shape parameter. The scheme predicts functions’ values at a set of dis- 
crete evaluation points, through a global sparse linear system. Compared to standard implementation, 
computational efficiency is achieved by using parallel computing. Applications of the proposed algorithms 
to 2D and 3D benchmark functions on uniformly distributed random points, the Halton quasi-points on 
regular or Stanford bunny shape domains, and an image interpolation problem confirmed the effectiveness 
of the method. We also compared the algorithms with other methods available in the literature to show 
the superiority of using PS augmented with a polynomial basis. High accuracy can be easily achieved by 
increasing the order of polyharmonic splines or the number of points in local domains, when small order 
of polynomials are used in the basis. MATLAB code for the 3D bunny example is shared on MATLAB 
Central File Exchange [1]. 
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1. Introduction 


Scattered data interpolation is a common and fundamental problem in many scientific and engineering 
studies [2-8]. There has already been so much work [9] in this area yet interpolation is still a difficult 
and computationally expensive problem. Polynomial interpolation and piece-wise polynomial splines have 
been generally used until 1971, Rolland Hardy introduced an interpolation method based on multiquadric 
(MQ) radial basis function (RBF) [10]. Since then, many different RBF interpolation schemes have been 
developed. A few different methods can be found in [11-13], and a comparison of radial basis functions 
methods can be found in [14]. 

Despite the simplicity, applicability to various kinds of problems, and effectiveness of RBF-based meth- 
ods, the resultant system of equations is often ill-conditioned when there is a large number of data points. 
Apart from that, global interpolation methods based on RBFs suffer from typical drawbacks of a global 
method such as high memory requirement and high computational cost. In fact, global RBF methods 
produce a linear system of size large as the number of data points. Generally, the costs of direct solution of 
such systems are O(N?) and O(N?) memory usage. There are several methods that get around these issues 
such as domain decomposition [15, 16], accelerated iterated approximate moving least squares [17], RBF- 
QR algorithms [18], a compactly supported RBFs [19], radial basis function-finite difference method [20], 
and many others. One disadvantage of the domain decomposition methods is that the domain discretiza- 
tion and joining phase of the local interpolants are not easy to implement. Recently, Hansen shared a 
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toolbox “regtools” including regtoolsTSVD, Tikhonov, and LSQR regularization methods for analysis 
and solution of discrete ill-posed problems [21]. 

The two methods proposed in [22] are fast localized RBF algorithms for large-scale 2D scattered data 
interpolation. In these methods, an interpolation is performed on each local influence domain and then 
all the influence domains are combined into a sparse matrix with the use of RBFs like Gaussian, MQ, 
Normalized MQ, and Matern. These methods are categorized as implicit methods in such a way that the 
interpolant is not explicitly defined by the numerical approximation, instead the approximation at sets 
of discrete evaluation points is given. We will continue to use such characterization throughout our work. 
On the other hand, if the interpolation function is defined by the numerical schemes, we call it an explicit 
method. 

The authors in [22] claimed that the conditionally positive definite RBFs such as polyharmonic splines 
r?°? ]nr and r2"+! “failed to produce acceptable results”. However, in this paper, we verified that when 
polyharmonic splines together with a polynomial basis is used on the implicit interpolation scheme, high- 
order accuracy, and excellent conditioning of the global sparse systems are gained without the need for 
selecting a shape parameter. The performance of the proposed algorithm is even more accurate than the 
most recent publication [8] in 2023. Although [8] is a global interpolation scheme, our algorithms are local. 

One may argue the proposed algorithm is the same as the idea of the Radial Basis Function — Finite Dif- 
ference Method (RBF-FD) [20] when applied to interpolation problems. However, there are fundamental 
differences between the two methods: 1) Our method is an implicit method, where only approximations 
at a set of discrete evaluation points are produced, but RBF-FD is an explicit interpolation in which 
an interpolation function is given by the numerical scheme; 2) Our method constructs the local domains 
by searching within the interpolation points or union of interpolation and evaluation points, but the 
RBF-FD constructs local domains purely within the evaluation points. 

In section 2, we briefly introduced global RBF interpolation and the positive-definiteness of the RBFs. 
In section 3, we propose the localized RBF interpolation methods for scattered data interpolation in 
Rĉ, where d is a positive integer. We chategorize this algorithm as a Local Implicit Interpolation using 
Polyhamonic Splines and Polynomials (LI2Poly2) Algorithm 1 and Algorithm 2 based on two different 
ways of creating local domains. Section 4 dedicates to numerical experiments carried on with 2D scat- 
tered data followed by a few 3D experiements. Numerical results are compared with the implicit local 
RBF method [22] and the CS-RBF method. The performance of the proposed method is demonstrated 
regarding accuracy, efficiency, and parameter selection. In section 5, we draw some conclusions on the 
usage of polyharmonic spline basis in the proposed method and possible improvements. 


2. Radial basis function interpolation 


The problem of scattered data interpolation is, given a set of distinct data points x; € R? with associated 
data values y; € R for i = 1,2,..., N, find an interpolant function f(x) : Rf > R, satisfying f(x;) = 
yi, 7=1,...,N. "i 

The global radial basis function interpolant f is given by 


N 
F(x) = $ aell- x;ll), (1) 
j=l 


where ¢(r) is a radial basis function with the r = ||x — x;|| > 0 defined as the Euclidean distance, and 
Xj is the center of the RBF. Note that r is the distance between point x and the centers of the basis 
functions. 

The unknown real coefficients a;,7 = 1,2,...,.N are determined by enforcing the interpolation condi- 


tion f(x:) = yi: 


N 
yi = flxi) = Yo oolli- xl), $= 1,2... (2) 
j=l 
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The resulting N x N linear system of equations can be represented by a matrix form 
Aa = b, 


where a = [a,Q2,...,an]", b = [y1, y2,-.-,¥n]", and entries of A are given by ai; = (||xi—x;||)), 1 < 
ij SN. 

The solution of the above interpolation problem exists and is unique, if and only if the interpolation 
matrix A is nonsingular, which is true for certain choices of RBFs that are positive definite. 


Theorem 2.1. [13] A real-valued continuous function @ is positive definite on R? if and only if it is even 
and 


N N 
dd, aiajọ(xi — xj) > 0, (3) 


for any N distinct data points x1,..., xy € R? and a = (@1,... ay)! € RX. 


Radially symmetric RBFs are clearly even functions. If ¢(r) is a positive definite function, it can be 
proved that the interpolation matrix A is a positive definite matrix for any distinct points X1,..., XN 
making it nonsingular. For example, Gaussian, inverse multiquadrics (IMQ), Matern, and compactly- 
supported RBFs (CS-RBFs) are positive definite functions [23, 24]. Other commonly used RBFs such 
as multiquadrics (MQ) and polyharmonic splines (PS) on the other hand have only been shown to be 
conditionally positive definite. 


Definition 2.2. [23] A real-valued continuous even function ¢ is called conditionally positive definite of 
order m on Rt if 


Sadik i— xj) > 0, (4) 


i=1 j=1 


ae . d = T N oh 
ee = weeds 
for any N distinct data points x, xy € R° and a = (a1 ay)’ ER” satisfying 


N 
S apx) = 0 (5) 
j=l 


for any real valued polynomial p of degree at most m — 1. The function @ is called strictly conditionally 
positive definite of order m on R? if the points x1,...,xn € R? are distinct, and a £ O implies strict 
inequality in (4). 


When using radial basis functions that are conditionally positive definite, one has to add polynomial 
basis functions of a certain maximal degree to the interpolate function (1) as follows 


N q 
x) =) aj 6(|x — xl) + $ Supe (x) (6) 
j=l k=1 


where {p1,p2,..-,Pq} forms a basis for Pe a space of polynomials of total degree less than or equal to 
m-— 1 in d-variables, where q = (Hm=, See Table 1 for some examples of the values of q in d dimension. 

The new linear system is created by enforcing the interpolation condition f G6) Sy Tort = NE 
We also consider the following additional insolvency constraints for the polynomial part to guarantee a 


unique solution for the new linear system: 


N 
Saype(xj) =0, k=1,..,4. (7) 


j=1 
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m-—1 | d=1 |d=2]|d=3 

0 1 1 1 

1 2 3 4 

2 3 6 10 
3 4 10 20 
4 5 15 35 
5 6 21 56 
6 7 28 84 

Table 1. The dimension of the basis {p1, p2,...,pq} where p;, i =1,...,q are polynomials of degree up to m — 1 in d dimension. 


Thus, the following linear system can be achieved by combining the interpolation condition and the 
additional insolvency constraints (7): 


T 
Coca! » 
P. O B 0 
where a = [a1,02,...,an]", B = [(1, Ba,-.. alts b = ly1,yo,---,yw]", and entries of A are given by 
aij = O(||xi — x;l|)), 1 < i,j < N and entries of P are given by pij = pj(xi), 1 <i < N1<j<¢. 
Once the unknown coefficients are obtained by solving the resulting (N + q) x (N + q) linear system, 
the interpolations can be made. Global interpolation methods based on RBF are easy to implement 
but, when the number of collocation points within the domain is high, the resulting matrix suffers 
from ill-conditioning and problems associated with computational cost including storage and time. On 
the other hand, adding polynomials in a global context can lead to other drawbacks associated with 
polynomial interpolation. With the limitations of global methods, formulations of localized methods offer 
an alternative for large-scale realistic data. In the next section, we present two local methods which have 
similar localization procedures as the generalized finite difference method. 


3. Localized RBF interpolation based on Polyharmonic Splines 


The two methods proposed in [22] are fast localized RBF algorithms for large-scale 2D scattered data 
interpolation. In these methods, an interpolation is performed on each local influence domain and then 
all the influence domains are combined into a sparse matrix with the use of RBFs like Gaussian, MQ, 
Normalized MQ, and Matern. All of these RBFs come with a free shape parameter that needs to be 
chosen carefully in the interpolation process. Yet, there is no practical and theoretical procedure for 
choosing the optimal shape parameter in applications except for some efforts made [25-28]. It is true that 
smooth RBFs like MQ give more accurate results at smaller values of the shape parameter c. But there 
is a trade-off between accuracy and conditioning. As c decreases function becomes flatter and accuracy 
increases until numerical ill-conditioning steps in. There are ways to improve the performance of MQ 
RBF methods, such as employing fictitious nodes [29, 30], pre-conditioning [16], etc [31-33]. 

We have improved two local methods in [22] by incorporating shape parameter-free polyharmonic 
splines (PS RBF) together with polynomial basis functions. This is an extension to the current methods, 
achieving high accuracy without the need of selecting a shape parameter. It is known that when high- 
order polynomials are used in global methods can lead to Runge’s phenomenon. Such effects are alleviated 
simply in the local context as one is only interested in interpolation within a small neighborhood. 

For convenience, let’s consider the interpolation problem in two-dimensional spaces. Suppose we have a 
set of distinct scattered data points {x;}®_, C R? and their function values { f(x;)}_, C R. Let T Cc 


R? be a set of evaluation points. We try to find the interpolant f such that f (z;) x f(z), j=1,..., Mi 
and f(x) = f(x), i= 1, see N. 
The polyharmonic splines in R@ are defined as follows: 


r?k-4]log(r), for d even 
þa k(r) = 2k-d 
r ‘ for d odd 


a where 2k > d, and with k € N. For example, in R?, ¢2.3(r) = r+ log(r), then d = 2, k = 3, thus polynomials 
4 up to degree 2 or higher need to be added to ensure the unique solvability. 


4 3.1. LI2Poly2—Algorithm 1 


so For each x;, we choose n nearest evaluation points to x; to create a local influence domain Q; = {zl l yes in 


sı which j = 1...n denotes for local indexing for each node in the Q;. Note that n < N. In this construction, 
52 each interpolation point has a local influence domain that contains only n evaluation points. Figure 1 
s shows an example of local domain construction with n = 5. 


@ Interpolation Points 
© Evaluation Points 


Figure 1. Local influence domain of x; with five nearest evaluation points in Algorithm 1. 


ti] 


7 = Q; for some j < n. 


Now let’s focus on the RBF interpolation on the local domain Q;, and let z = z 


Thus, f(z) can be written as previously shown in (6) which is the following, 
n ; q 
fœ) = X oeli - 251) + J ony pil). 
j 
j=1 l=1 


154 where ¢ is chosen to be polyharmonic splines RBF and {pj,p2,...,Pq} forms a basis for polynomials up 
155 to degree less than or equal to m — 1 in R, q = eae ee Notice that m = k — [d/2] + 1 is the order of 
156 the conditionally positive definite polyharmonic splines ġa, [34]. In addition, the size of the local domain 
1s7 should be greater than the number of polynomial basis functions q = Ge) for this to work. In this 
iss example, the local domain size needs to be larger than q = CE) = 6 in R?. 


2 
Collocation on the local domain of influence leads to the following system: 


Y a; (Ila! — z?) + So ong vila!) = fiat), &=1,2,...n, (9) 
j=l {=l 
Y oap) =0, 1=1,2,...4. (10) 
j=1 


159 Let the coefficient matrices on the first and second terms of the left-hand side of (9) be ® and P 
160 respectively. That is 


-1 
i ld l) of l zi y Yi 
m b(||Zy° — zi ||) (lz — zn ||) and P= a1) 
1 ta Yn yo 
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Then we can introduce a block matrix from for the system (9)-(10) as follows: 


(E E)oH-(4) e 


where al = [ay,ae,...,Qntq]", f, = LÊ), F(ab), Sui fay. Let the coefficient matrix in (12) be 
W. Then the unknown coefficients in (9)-(10) can be expressed as 


all = y7! ( - i (13) 


Therefore, for i = 1,..., N, 


n , q 
Fœ = X ajo (|| — 2l) + Y ant pili), 
j=1 l=1 


SN Ga fo (14) 
where 


(x;) = bie =A easke- a e a g N 


and An+g(xi) = ®(x;)W'. Note that A,(x;) is obtained by omitting last q elements of the vector 
Antq(Xi). re 

The system (14) can be reformulated to an N x N; sparse system easily by extending local fn to a 
global fy, = [f(z1), f(z2),...,f(zn,)|". This can be done by inserting zeros to An(x;) based on the 
mapping between f, and În,- It follows that 


a 


f (xj) = An, (xi) fn, (15) 


where Ay, (x;) is a vector that is obtained by adding zeros into A,,(x;) appropriate at places. 

For example, assume N; = 50,n = 3, and Q; = {all oll ofl = {Z12,Z20, 223}. Then we need to 
insert 47 zeros into An(x;), while only keeping non zero values at 12%, 20¢ and 23”? positions. i.e, 
An, (xi) = [0,0,..., # ,0,0,...,0, # ,0,0, # ,0,0,...,. 0 ]. By solving (15), the unknown approx- 

Ww Ww Ww aos 

12th 20t} 23rd 50%% 
imation function values at N; evaluation points, f (zj), j = 1,2,..., N; can be found. This can be done 
using any direct or iterative methods for solving systems of linear equations given that N < N. The 
proposed algorithm works for any dimension d like many other RBF-based methods. 


3.2. LI2Poly2—-Algorithm 2 


Algorithm 2 is developed with the idea of taking all the interpolation points and evaluation points in 
each local neighborhood into consideration for the interpolation at a single point. In this construction, 
each interpolation point will have a local influence domain that contains both evaluation points and 
interpolation points, total into n. Figure 2 shows an example of local domain construction with n = 7. 
For each x;, we choose the nearest nı evaluation points and nz interpolation points which adds up to n 


to create a local influence domain 2; = piri Dea 1- The local interpolation procedure presented 


183 


184 


185 


186 


187 


188 


189 


190 


0. 9 i 
@0e © © o O 
O ö O 
oa eee @ e- 
0.7 s O O 
© eO è o e 
0.6 (a) 
e oOo o @ 
oF O o 
© © o e 
0.4 
(e 
E & @ © e 
; k m 
Moe © 000 e 
© O 
0.1 "a e o e %® Ce (J 


0 of o2? 03 04 OF Of 07 NR Na 4 
@ Interpolation Points 
O Evaluation Points 


Figure 2. Local influence domain consists of seven nearest points to x; in Algorithm 2. 


in Algorithm 1 resulted in (15) should be changed to account for the changes in Algorithm 2. The resulting 
system of linear equations is size N x (N; + N) and is underdetermined as follows, 


a 


Hsin ( > ) de (16) 


where fn = [f(xi),...,f(xw)]. In order to solve this system, the system has been reformulated to a 
2N x (Ni + N) system as follows 


Onxn, Inxn fN fN j` 
The approximated function values fn, can be obtained by solving (17) using any system-solving technique. 


kd-tree amoung Nt evaluation points O(Nilog(N:)) 
Local domains of N interpolation points O(nNlog(N:)) 


Solve N small local systems of (n+q)x(n+q) O(N(n + q)°) 


Solve a sparse systemof NxN O(Nn?) 


kd-tree amoung Nt+N points O((N + N:)log(N + Nz)) 


Algorithm 2 
O(nNlog(N + N:)) 


Local domains of N interpolation points 


O(N(n4 q)’) 


Solve N small local systems of (n+q)x(n+q) 


Solve a sparse system of (N+Nt)x(N+Nt) © (OY + Nv) n?) 


Figure 3. Computational complexity in time with regards LI2Poly2 Algorithm 1 and Algorithm 2. 


The efficiency of the method can be illustrated by considering its asymptotic computational complexity. 
For both algorithms, we need to 


Step 1. For each interpolation point, we need to calculate the coefficient matrix in (14) by 
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o first, build kd-tree amung the all N; evaluation points for Algorithm 1 or for all N; evaluation 
points and N interpolation points for Algorithm 2; 


o second, find n nearest neighbors of all N interpolation points using the kd-tree created above; 
o third, solve small local systems of size (n + q) x (n + q), where there are N such systems; 
Step 2. Solve a sparse system of size N x N, in which n unknowns in each equation. 


Figure 3 shows detailed computational complexity associated with both algorithms. Note that there 
are N, evaluation points, N interpolation points, and q is the number of polynomial basis. Thus, there 
are N small systems of size (n+ q) x (n + q) in Algorithm 1. The cost for solving these N small systems 
is O(N(n + q)?). However, n,q << N in practice. Hence, the computational complexity in time for 
Algorithm 1 is O((N; + nN)log(Nz)) + O(N (n + q)?) + O(Nn?) and the computational complexity of 
Algorithm 2 is O(((n + 1)N + Nz)log(N + Nz)) + O(N (n + q)?) + O((N + Nt)n?). 

In the next section, we demonstrate the accuracy and efficiency of the proposed method using two- 
dimensional data intensively and some results from the three-dimensional data as well. The stability of 
the method is also inspected by studying the condition number of the interpolation matrices and global 
sparse matrices. 


4. Numerical Results 


To illustrate the effectiveness of the method on scattered data, we consider examples of two and three 
dimensions on both regular and irregular domains. Recall the following key parameters/notations: 


N : the number of interpolation points 

N;: number of evaluation (test) points 

n : the number of points in the local domain of influence 
m : the degree of highest-order polynomials 

k : the order of PS. 


Numerical results are compared with the exact function values and the accuracy of the method is 
measured in terms of root mean squared error (€rms) and the maximum absolute error (€maz), whose are 
given by 


N 
1 P N 3 
Erms = N X (fi = fd, Emax = max] fi E fil (18) 
i=l 


where f; = f (xi) is the approximated value of f; = f(x;) . The condition number of the interpolation 
matrix A is defined as 


cond(A) = ||All2||4~I|2 = 7 (19) 


where Cmax and Omin are the largest and smallest singular values of A. In our numerical experiments, 
we choose the interpolation points to be evenly distributed, while the evaluation points (test points) 
are randomly distributed Halton quasi-points with the constraint of N; < N. In a case of N < Ni, we 
split the evaluation points into subsets and perform the interpolation separately. We will explain how 
interpolation is done in such cases later in this section. 

All numerical experiments have been performed in MATLAB on a MacBook Pro with a 3.2GHz Apple 
M1 processor and 16 GB memory. The algorithm code has been parallelized using the Matlab Parallel 
Computing Toolbox for improving the performance. Moreover, construction of the local domains by 
searching the nearest evaluation points has been done using the Matlab built-in function knnsearch from 
the statistics and machine learning toolbox. 


Example 4.1. In the first example, we investigate the performances of the proposed methods on Franke’s 
six test functions [85] on the unit square [0,1] x [0,1]. For simplicity, we mainly show details on the test 
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function F,, the algorithms behave in a similar way on all other five test functions. The test function Fi 
is provided as follows. 


$ exp ( } ((9a 
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a —Alg1 
10°19; a =-=- Alg 2 | 
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(b) Rate of convergence w.r.t spatial variables. 


(a) Profile of test function F3; 


Figure 4. Franke’s benchmark test function F; [35] on the left, and the root mean squared errors, €rms, versus average distance 
between nodes, h, for Fi. 


The left of Figure 4 shows the profile of test function Ff. On the right of Figure 4, the rate of convergence 
is displayed for the two proposed algorithms and the reference implicit method, with respect to spatial 
discretization. From the figure, it can be seen that Algorithm 2 has the highest convergence rate with 
the highest accuracy. 


One of the main advantage of PS is that it does not have a shape parameter that needs to be determined 
during the interpolation process. However, we still need to select other parameters such as the order of 
the PS RBF (k), the order of the polynomial basis (m), and the number of local points (n). Since 
both algorithms behave similar and Algorithm 2 performs better in terms of accuracy without loss of 
computational efficiency, we examine the performance of Algorithm 1 (the worse case senarior) as we 
vary k,m, and n. Note that the findings presented here are comparable to those of Algorithm 2. 


7 22 
10 5 r r r r 25 x10 : : : : 

—n=20 @ 
—n=30 8 

— n=40 z | 
10°F J 5 
5 
é 2 

o 3 7 
W 27, | Ē 
n 10 E 
= 2 

oc pool 4 
fo} 
fo} 
E. 
10°F 7 5 

= | 
= 
O 
oO 
x 
ð 
10° i = 

0 0.5 15 2 5 Bs 
N x10f x10 


Figure 5. RMS errors versus N (left) and maximum condition number of local interpolation matrices versus N (right) for various n 
on Fı with k = 4,m = 3 in Example 4.1 using Algorithm 1. 


From the left of Figure 5, we can see that interpolation error decreases as the number of interpolation 
points increases, and three separate lines indicate that one may obtain even more accurate results by 
increasing the number of points in the local domains. From the right of Figure 5, we observe that the 
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ill-conditioning grows as the number of interpolation points increases, although the error decreases with 
increasing N as shown in the left of Figure 5. 
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Figure 6. Left: maximum errors versus the order of PS k for N = 100?, Nz = 9000, n = 30 and m = 6 using different test functions 
and Algorithm 1. Middle: maximum errors versus the polynomial order m for N = 100? , N; = 9000, n = 30 and k = 4 using different 
test functions. Right: CPU time (in seconds) versus m for test function Fy using Algorithm 1 in Example4.1. 


The left of Figure6 indicates that the accuracy doesn’t change significantly when the order of PS 
changes (about one order of magnitude differences). In the middle of Figure6, it indicates that when 
the order of the polynomial basis m increases, the accuracy improves. Yet, one would not need to in- 
crease the order of the basis too much as it will lead to a high computational cost as evident from the 
right of Figure6. However, we can reduce the simulation time significantly by using parallel comput- 
ing while calculating the weights associated with each node in the global sparse system. For instance, 
when using k 4,m 6,N = 40,000 for Fi, the CPU time is reduced from 1050 seconds to 250 
seconds. Note that the results from Franke’s test functions F4 and F; [35] are also included, and the 


test functions are presented below for reference: Fy(x,y) = Lexp | SI (œ ie | (y n’) and 


4 2 2 
comparison with F}, the accuracy of interpolation results from those two functions is much better than 


those from F. 


F(x,y) = = exp si (@ 1j” + (y — p^ . Due to the relative smoothness of those two functions in 
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Figure 7. Sensitivity respect to points distribution: RMS errors versus p with N = 100?, N; = 90?, and n = 30 using PS of order 
k = 4 and polynomial basis of order m = 6 in Example 4.1. 


To validate the stability of the methods, we examine the sensitivity of the algorithms with respect to 
the locations of the points. Here we consider the perturbed grid points provided by 
fi = zi + paz Ôx 


yi + py" dy 


Yi 


where (x; yi) € (0,1)? are uniformly distributed grid points. Let p denote the degree of randomness 
and 6x, ôy are the shortest distance between two points in x and y directions. Lastly, ee aye and are 


uniformly distributed random numbers in (0, 1). Figure 7 shows the RMS error obtained using Algorithm 1 


10 


and Algorithm 2 on the left and on the right, respectively, for various degrees of randomness. When p is 
zero, the points are evenly distributed and as p increases, the interpolation points become more distorted. 
However, the error does not change much for all three test functions. Hence, the stability with regard to 
the point distribution is validated. 

Table 2 shows the condition number of the global sparse matrix, and the maximum condition number 
of all local collocation matrices when using PS of order k = 4 and m = 6, n = 30 on F} for two 
algorithms. This case was chosen mainly because this is a general case when reasonable accuracy and 
efficiency are achieved. Clearly, the small collocation matrices have higher condition numbers when the 
number of interpolation points and evaluation points increases, even becoming ill-conditioned in many 
cases. However, the global sparse matrices always have excellent condition numbers regardless of which 
algorithm and what parameters are used. When dealing with ill-conditioned small local systems, we would 
recommend a singular value decomposition, TSVD, or pre-conditioning techniques [21]. 

Furthermore, when comparing the CPU time, Algorithm 1 demonstrates greater efficiency compared to 
Algorithm 2. By examining the CPU time of the proposed methods, the most time-consuming part of the 
algorithms is identified as the construction of the local matrices while the final sparse system solving time 
is negligible. The algorithm knnsearch in MATLAB for nearest-neighbor classification is not the most 
efficient such algorithm. The highly efficient kd-tree algorithm [36] can be used if desired computational 
efficiency is needed. 


| Algorithm 1 Algorithm 2 

| (N, Ni) (502,2000) | (1507, 20000) (502,2000) | (1507, 20000) 
| Condition No. of sparse matrix 2.89E+03 4.81E+04 3.78E+03 8.13E+06 

| Maximum condition No. of local matrices | 1.89E+14 2.8405E+20 | 3.96E+18 1.04E+22 

| CPU time (s) 1.71 75.60 3.09 166.37 


Table 2. Condition numbers and CPU time for Algorithm 1 and Algorithm 2 using PS of order k = 4 and m = 6, n = 30 on F; in 
Example 4.1. 


The algorithms introduced in this paper are very much similar to the global RBF interpolation methods 
using CS-RBF. The Wendland’s CS-RBFs were first introduced in [37], which is piecewise polynomial basis 
such as ¢(r) = (1 — r)4.(4r + 1). Interpolation with CS-RBFs globally also resulted in a global sparse 
system [38]. Thus, we examine our algorithms against CS-RBFs interpolation. It is worth mentioning 
that the CS-RBF interpolation is still a global method since the interpolation matrix is created using 
CS-RBF's with all interpolation points as the centers. Table3 shows performance in terms of accuracy 
and efficiency. In this example, PS of order 4 is used with a polynomial basis of order 3. From previous 
analysis, we know that we can improve the accuracy using a higher-order of polynomial basis or more 
neighboring points. With slightly higher accuracy, we still able to be much more efficient in terms of 
computational time. 


N N CS-RBF Proposed Algorithm1 | Proposed Algorithm2 
t Ermis CPU time Ermia CPU time Erms CPU time 
1007 | 9000 | 7.66E—07 186.0 6.00E—08 15.10 4.02E—08 30.27 
150? | 20000 | 2.01E—07 1454.8 8.34E—09 71.2 5.26E—09 136.57 
Table 3. Comparison of CS-RBF method and the proposed methods with k = 4,m = 3,n = 30 for function Fı in Example 4.1. 


Table 4 shows that the interpolations achieved with PS RBF have better accuracy than the other RBFs, 


Algorithm 1 is utilized to interpolate the test function Fg = z [64 — 81 (e — Ja + (y — 3’) — 5. 


The Leave-One-Out-Cross-Validation method (LOOCV) [26] is employed to determine the best shape 
parameter associated with RBFs with shape parameters to ensure that we choose the “optimal” value to 
the best of our ability. In LOOCV, we aim to fit data several times using a different training set (all but 
one data point) and testing set (one of the data points) each time, then calculating the test root mean 
squared error to be the average of all of the test. 
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RBF Erms Emax C 
Gaussian 2.15E—13 | 4.06E—12 | 4.994 
MQ 1.17E—12 | 2.97E—11 | 3.208 
IMQ 3.29E—13 | 9.93E—12 | 2.629 
Matern(order 2) | 1.71E—12 | 3.43E—11 | 1.909 
PS4 3.22E—14 | 6.26E—13 - 


Table 4. Comparison of erms and €max using various RBF with N = 1002, N: = 9000 and n = 30 for Fg Example 4.1. 


Example 4.2. In this example, we further investigate the performance of the proposed Algorithm 1 using 
the function: 


300 


301 


FT(x,y) = Vx? +y? + 0.2, 


which has a singularity and is more challenging to interpolate. 


(20) 


302 


x104 


0 


Figure 8. Profile of test function Fy on the left and absolute errors of interpolation with N = 150?, N¢ = 20,000,n = 30, k = 4,m = 6 
on the right in Example 4.2. 


Figure 8 shows the profile of F in [—1,1] x [—1,1] on the left, and a profile of absolute error for test 
function Fy on the right, where k = 4, m = 6, N = 1502, N; = 20,000 and n = 30. The absolute error 
is reasonably accurate, and it can be improved by employing more points in the local domains. The 
proposed method can certainly be successfully used for the interpolation of challenging functions without 
implementing any special treatments such as grid refinements. 
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Figure 9. Set of evaluation points and consistently distributed two sets of interpolation points. 


308 We noticed that the higher the number of evaluation points, Ny, the greater the accuracy of the 
39 algorithms, but N < N. However, when there N; > N, we can extract sub-samples of dimension Ni << N 
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which can be interpolated group by group using the proposed method. These subsets need to be extracted 
as uniformly as possible from the domain to avoid any discontinuity between interpolated values. We 
implement this for the function Fy with N = 22,500, and various N; that are higher than the given 
N. As shown in Table 5, the proposed method can maintain a similar level of accuracy for all the cases 
despite the number of evaluation points. An increase in the CPU time is expected as the number of 
sub-samples gets higher. However, this algorithm is suitable for parallel implementation as these groups 
are chosen independently, thus computer running time can be reduced significantly. 


Nt No. of €rms CPU time 
subsets 
30,000 | 2 2.70E—05 126.32 
40,000 | 2 2.19E—05 151.40 
50,000 | 3 2.14E—05 196.44 
60,000 | 3 2.06E—05 220.11 


Table 5. RMS errors and CPU time of Fy interpolated with N = 22,500,n = 30 and various N+ (N < N+) in Example 4.2. 
Example 4.3. Interpolation is one of the common methods used to enhance image quality. In this ex- 
ample, 256 x 256 pixels 2D gray-scale Lena image shown on the left of the Figure 10 is enhanced. The 
enhanced 512 x 512 pixels image is shown on the right. Note that the enhanced image is obtained by using 
Algorithm 1 with PS of order k = 1, a polynomial of order m = 0, and n = 3. As the set of evaluation 
points is higher than the interpolation points in this problem, interpolation is performed group-wise. By 
maintaining a low order of polynomial basis and restricting the size of the local domains, we achieved 
visible improvements to the image while keeping the computational time low. 


50 100 150 200 250 300 “350 400 450 500 
Figure 10. Profile of the original Lena image on the left and the enhanced Lena image on the right in Example 4.3. 


Example 4.4. The purpose of this example is to test Algorithm 1 on three-dimensional space. We use a 
computational domain of the Stanford Bunny, which is more challenging and irregular. The boundary data 
points for this domain are available at the website of the Stanford Computer Graphics Laboratory [39]. 
Algorithm 1 was tested on function H in (21), which is a trivariate test function used in/14]. The function 
plotted on the bunny surface can be found in Figure 11 (right): 


H(z,y,z) = : ((2 0.5)? + (y — 0.5)? + (z— 0.5)*)] (21) 
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In our numerical simulations, N = 4234 of interpolation points which consist of 2345 interior points 
and 1889 boundary points are used to interpolate N; = 1143 of test points. The computational domain 
with interpolation points lying on the boundary surface is shown in Figure 11 (left). Since the original 
bunny data scale is too small for this application, the coordinates of bunny points are multiplied by 10. 
Algorithm 1 is employed with PS of order k = 4, the polynomial basis of order m = 3, and n = 55. Figure 
12 (left) shows the profile of the absolute errors on the surface of the Bunny and Figure 12 (right) further 
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Figure 11. Profile of the computational domain (Stanford Bunny) with the set of interpolation points on the boundary (left) and the 
profile of H function on the surface of the bunny (right) in Example 4.4. 


confirms that the absolute errors in the interior of the bunny remain as low as the order of 5. The largest 
absolute error in all evaluation points is within an order of 1075, with the majority of evaluation points 
having absolute errors less than 1076 (which can be found in the tallest bar in the histogram). It verifies 
the method’s high accuracy despite the complexity of the domain. 
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Figure 12. The profile of the absolute errors on the surface of the bunny (left) and frequency or count of data points in the interior 
of the bunny whose absolute error are within each range (right) in Example 4.4. 
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Figure 13. The maximum absolute error, €max versus parameter n on the left, and corresponding CPU time on the right in Exam- 
ple 4.4, where N = 4234, N; = 1143. 


On the left of Figure 13, it shows the maximum absolute errors, €maz, versus the number of points in 
local domains, n. Note that the range of the number of points in local domains is large, spanning from 25 
to 1000. Regardless of the size of n, the algorithm consistently improves in accuracy. Thus, the greater 
the number of points in local domains, the higher the accuracy of the algorithm. However, we recommend 
choosing n to be less than 200, as larger values result in increased computational time (as shown on the 
right of Figure 13). 


Figure 14. The maximum absolute error, €maz versus parameter m on the left, parameter k on the right in Example 4.4, where 
N = 4234, Nt = 1143. 


Figure 14 shows the maximum absolute errors using Algorithm 1 when 4234 interpolation points and 
1143 evaluation points are used in Example 4.4. The left of Figure 14 plots em, versus the polynomial 
order m when the number of points in local domains is chosen as n = 85 and n = 120. It can be seen 
from the figure that the order of the polynomials basis, m should not be too large. The algorithm easily 
achieves an accuracy of order 5 when m is as small as 2. 

On the right of Figure,14, the plot displays €m¢, versus the order of PS, k when the number of points 
in local domains is chosen as n = 25,45, and n = 85. Please note the following: (1) as we discussed 
before, increasing the number of local points improves the accuracy, hence we did not present results for 
n = 120; (2) m = 3, n = 45 is sufficient to achieve reasonable accuracy, so the effect of the order of PS 
is not evident; (3) when n = 25, k has to be small. In this case, the best accuracy is obtained only when 
k=l; 

Thus, we recommend using enough number of points in local domains to achieve high accuracy while 
keeping the order of polynomial basis small. In addition, the order of PS does not need to be high as it 
does not significantly affect the accuracy when a sufficient number of local points is used. 

MATLAB code for this example is shared on MATLAB Central File Exchange [1]. 
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5. Conclusion 


In this paper, we improved two implicit localized RBF interpolation methods using polyharmonic splines 
and a polynomial basis for scattered data interpolation that was proposed in [22]. The method uses 
evaluation points or a combination of evaluation and interpolation points as the search domain to create 
local domains of influences for each interpolation point. The resulting linear system is a sparse system 
with the evaluation points’ function values as the unknowns. The original paper claims the method “does 
not produce reasonable results when using the polyharmonic splines”. However, we discovered in this 
paper, the claims are not true. 

The interpolation with polyharmonic splines and low order polynomials can be solved with great 
accuracy and efficiency. The higher the order of polyharmonic splines or the number of points in local 
domains, the better the accuracy becomes, under the assumption that the points in local domains are 
sufficient. Detailed conclusions are drawn below: (1) RBFs such as polyharmonic splines and polynomials 
contain no shape parameter and produce higher accuracy compared to other RBFs; (2) Algorithm 1 is 
based on the construction of local influence domains entirely from the evaluation points and Algorithm 2 
takes both interpolation and evaluation points into consideration. Both algorithms were found to be very 
attractive, easy to use in 2D and 3D problems, and mostly able to overcome the downsides of global RBF 
interpolation using polyharmonic splines and polynomial basis; (3) The computational complexities of 
Algorithm 1 and Algorithm 2 are very close. Additionally, Algorithm 2 is more accurate than Algorithm 
1 but we mainly focused on Algorithm 1 due to it’s minor improvement of efficiency. 

In summary, our proposed interpolation algorithm can deal with high-dimensional data interpolation on 
complicated domains, without the hassle of searching for shape parameters and without loss of accuracy 
and can be parallelized easily. This is a great improvement in terms of computational efficiency. Future 
work concerns deeper theoretical analysis of the convergence analysis and error estimates. In addition, 
we aim to to further reduce the computational complexity of Algorithm 2 and apply it to classification 
problems. 
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